                                                      # -*- coding: utf-8 -*-
"""
Created on Tue Sep 29 02:36:52 2020

@author: sarupa
"""
from __future__ import division
import numpy as np

# #==========================================#

Fp=0.5
T10=300
T20=300
V1= 4.0
V2= 4.0
V3= 4.0#4.0

E1= 50000
E2= 60000
k1= 9.972e06
k2= 9.36e06

dH1= -6*10000
dH2= -7*10000
aA= 3.5
aB= 1
aC= 0.5
Cp= 4.2
R= 8.314
rho= 1000
xA10=1.0
xA20=1.0
xB10= 0.0
xB20= 0.0
Hvap1= -3.57*10000
Hvap2= -1.57*10000
Hvap3= -4.07*10000
cm= 2


Nx=9
Nu=6
Ny=9
# Create CSTR ode model
#==========================================#
def cstr_xy_ode_scale(x_scale, u, w):
   
    xA1 = x_scale[0]
    xB1 = x_scale[1]
    T1 = x_scale[2]
    xA2 = x_scale[3]
    xB2 = x_scale[4]
    T2 = x_scale[5]
    xA3 = x_scale[6]
    xB3 = x_scale[7]
    T3 = x_scale[8]
        
    xC3 = 1 - xA3 - xB3
    x3a = aA*xA3 + aB*xB3 + aC*xC3
    xAr = aA*xA3/x3a
    xBr = aB*xB3/x3a
    xCr = aC*xC3/x3a
    
    
    
    F10= u[0]
    F20= u[1]
    Q1= u[2]
    Fr=u[3]
    Q3= u[4]
    Q2= u[5]
    
    F1 = F10 
    F2 =F1 + F20
    
    dx_scaledt= [ F10*(xA10 - xA1)/V1 - k1*np.exp(-E1/(R*T1))*xA1 +w[0],
         F10*(xB10 - xB1)/V1 + k1*np.exp(-E1/(R*T1))*xA1 - k2*np.exp(-E2/(R*T1))*xB1 +w[1],
         F10*(T10 - T1)/V1- dH1*k1*np.exp(-E1/(R*T1))*xA1/Cp*cm/rho - dH2*k2*np.exp(-E2/(R*T1))*xB1/Cp*cm/rho + Q1/(rho*Cp*V1) +w[2],
         F1*(xA1 - xA2)/V2 + F20*(xA20 - xA2)/V2 - k1*np.exp(-E1/(R*T2))*xA2 +w[3],
         F1*(xB1 - xB2)/V2 + F20*(xB20 - xB2)/V2 + k1*np.exp(-E1/(R*T2))*xA2 - k2*np.exp(-E2/(R*T2))*xB2 +w[4],
         F1*(T1 - T2)/V2 + F20*(T20-T2)/V2 - dH1*k1*np.exp(-E1/(R*T2))*xA2/Cp *cm/rho - dH2*k2*np.exp(-E2/(R*T2))*xB2/Cp *cm/rho+ Q2/(rho*Cp*V2) +w[5], 
         F2*(xA2 - xA3)/V3 - (Fr+Fp)*(xAr - xA3)/V3 +w[6],
         F2*(xB2 - xB3)/V3 - (Fr+Fp)*(xBr - xB3)/V3 +w[7],
         F2*(T2 - T3)/V3 + Q3/(rho*Cp*V3)+(Fr+Fp)*(xAr*Hvap1+xBr*Hvap2+xCr*Hvap3)/(rho*Cp*V3)*cm +w[8]
         ]
    return dx_scaledt
#==========================================#
# measurement functions 
def measurement_xy_model(x):
    C_matrix= np.zeros((Ny,Nx))

    C_matrix = np.array([[0.,0.,1.,0.,0.,0.,0.,0.,0],
                         [0.,0.,0.,0.,0.,1.,0.,0.,0],
                         [0.,0.,0.,0.,0.,0.,0.,0.,1.]]) 
    
    y_3 = np.dot(C_matrix,x)
    return y_3


